#include "stdafx.h"

//*****************************************************************************
int initTemperatureField(material *substance, int **grid, int row, int col, 
						 double tempEdge, 
						 double **tempCurr, double **tempNext) 
{
	int k;
	for (int i = 0; i < row; i++)
		for (int j = 0; j < col; j++)
		{
			k = grid[i][j];
			//if (k == EDGE_ID)
			//	tempCurr[i][j] = tempNext[i][j] = tempEdge;
			//else
				tempNext[i][j] = tempCurr[i][j] = substance[k].initTemp;
		}
	tempCurr[0][0] = tempNext[0][0] = substance[0].initTemp;

	return 1;
}

//*****************************************************************************
double calcTimeStep(material *substance, int **grid, int row, int col, 
					double htcInter[NUM_MATERIAL][NUM_MATERIAL], double dx) 
{
	double timeMax = 100000, timeTemp;
	double dxdx = dx * dx;
	for(int i = 1; i < row - 1; i++)
		for(int j = 1; j < col - 1; j++)
			//if (grid[i][j] != EDGE_ID)
			{
				timeTemp = 
					(substance[grid[i][j]].density 
					* substance[grid[i][j]].sh 
					* dxdx)	
					/ (htcInter[grid[i][j]][grid[i - 1][j]] 
					+ htcInter[grid[i][j]][grid[i + 1][j]] 
					+ htcInter[grid[i][j]][grid[i][j - 1]] 
					+ htcInter[grid[i][j]][grid[i][j + 1]]);
				if (timeTemp < timeMax)
					timeMax = timeTemp;
			}
		return timeMax;
}

double fixTimeStep(double maxTimeStep, double factor) 
{
	return maxTimeStep * factor;
}

//*****************************************************************************
int calTemperationNextStep(material *substance, int **grid, 
						   double htcInter[NUM_MATERIAL][NUM_MATERIAL],
						   int i, int j, double timeStep, double dxdx, 
						   double **tempCurr, double **tempNext) 
{
	int k = grid[i][j];
	tempNext[i][j] = tempCurr[i][j] + 
		timeStep / substance[k].density / substance[k].sh / dxdx * 
		(htcInter[k][grid[i][j - 1]]*(tempCurr[i][j - 1] - tempCurr[i][j])+
		htcInter[k][grid[i][j + 1]]*(tempCurr[i][j + 1] - tempCurr[i][j])+
		htcInter[k][grid[i - 1][j]]*(tempCurr[i - 1][j] - tempCurr[i][j])+
		htcInter[k][grid[i + 1][j]]*(tempCurr[i + 1][j] - tempCurr[i][j]));
	return 0;	
}

double fourierIteration(material *substance, int **grid, int row, int col, 
						double htcInter[NUM_MATERIAL][NUM_MATERIAL], 
						double timeStep, double dx, terminationCondition2D tc, 
						double **tempCurr, double **tempNext) 
{
	double	**swap;
	double	dxdx = dx * dx;
	double	cycle = 0;
	outputTemperatureField(tempCurr, row, col, cycle);
	printf("center of plate - %f @ time - %f\n", tempCurr[tc.y][tc.x], cycle);
	while(tempCurr[tc.y][tc.x] > tc.terminationTemp)
	{
		cycle += timeStep;
		for(int i = 1; i < row - 1; i++)
			for(int j = 1; j < col -1; j++)
			{
				calTemperationNextStep(substance, grid, htcInter, 
										i, j, timeStep, dxdx, 
										tempCurr, tempNext);
				//tempCurr[i][j] = tempNext[i][j];
			}
		swap = tempCurr; tempCurr = tempNext; tempNext = swap;
		if ((int)(cycle / OUTPUT_FREQUENCY) > (int)((cycle - timeStep) / OUTPUT_FREQUENCY)) 
		{
			outputTemperatureField(tempCurr, row, col, cycle);
			printf("center of plate - %f @ time - %f\n", tempCurr[tc.y][tc.x], cycle);
		}
	}
	outputTemperatureField(tempCurr, row, col, cycle + OUTPUT_FREQUENCY * 2);
	printf("center of plate - %f @ time - %f\n", tempCurr[tc.y][tc.x], cycle);
	return cycle;
	
}

//*****************************************************************************
//fourier method realization
double fourier2D(material *substance, int **grid, int row, int col,
				 double htcInter[NUM_MATERIAL][NUM_MATERIAL], 
				 double tempEdge,  double dx, terminationCondition2D tc)
{
	double		**tempCurr, **tempNext;
	tempCurr = (double **)malloc(sizeof(double *) * row);
	tempNext = (double **)malloc(sizeof(double *) * row);
	for(int i = 0; i < row; i++)
	{
		tempCurr[i] = (double *)malloc(sizeof(double) * col);
		tempNext[i] = (double *)malloc(sizeof(double) * col);
	}
	//timer
	LARGE_INTEGER	nStartCounter, nStopCounter;
	QueryPerformanceCounter(&nStartCounter);//timer starts

	initTemperatureField(substance, grid, row, col, tempEdge, tempCurr, tempNext);
	//outputTemperatureField(tempCurr, row, col, 0);

	QueryPerformanceCounter(&nStopCounter);//timer ends
	double nTime1 = getCostTime(nStartCounter, nStopCounter);

	double maxTimeStep = calcTimeStep(substance, grid, row, col, htcInter, dx);
	printf("max time step = %f \n", maxTimeStep);
	double timeStep = fixTimeStep(maxTimeStep, TIME_STEP_FACTOR);
	printf("time step = %f \n", timeStep);

	QueryPerformanceCounter(&nStopCounter);//timer ends
	double nTime2 = getCostTime(nStartCounter, nStopCounter);

	double cycle = fourierIteration(substance, grid, row, col, 
									htcInter, timeStep, dx, tc,									
									tempCurr, tempNext);

	QueryPerformanceCounter(&nStopCounter);//timer ends
	double nTime3 = getCostTime(nStartCounter, nStopCounter);

	return cycle;
}